1.1. Surface Matching : example on human scapula¶

This notebook presents a method to match a point cloud taken with any scanning system to a digital 3D scan in stl format.

1.1.1. Introduction¶

1.1.1.1. Method high lines¶

Surface matching is based on an optimization strategy. The objective of this optimization is to identify the parameters of the reference change allowing to pass from the coordinate system of the STL file to the coordinate system of the scan system.

This change of reference frame is formalized by a passage matrix formed by a rotation vector and a translation vector according to the Rodrigues formalism. There are thus a total of 6 DOF.

The optimization parameters of the cost function will therefore be these two vectors. The residual is defined by the signed distance between the points coming from the measurement systems after the change of reference frame and the stl mesh. This distance is calculated between the points of the closest scanning system and the STL mesh.

This problem being badly formulated, it is necessary to establish an initial solution in order to have a correct result.

This is why in a first step, we establish an initial using points that are identified on the STL model but also on the real model. A minimum of 3 points are necessary to block the 6DOF.

Then, in a second step, an optimization will be performed with the whole scan taking as initial solution the one previously calculated.

1.1.1.2. Modules importations¶

import trimesh
import numpy as np
from scipy import optimize
import plot_mesh as pltm
import plotly.graph_objects as go
import pandas as pd
import gbu

1.1.1.3. Case of study¶

Left Human scapula stl file

fig = pltm.create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")
fig.show()

1.1.2. Define reference points on STL file¶

We define reference points on the STl that can be easily found on the real model A point on the glenoid, coracoid and acromion is selected.

p_glene_stl = np.array([0.02515289854664963, 0.03122873596800735, 0.036798809060995745])
p_coracoide_stl = np.array(
    [0.026787610816186087, 0.05609095153197788, -0.0043762069034760384]
)
p_acromion_stl = np.array(
    [-0.008420164259823573, 0.07140986464971243, 0.018059532550748165]
)
fig = pltm.add_points(points=p_glene_stl, fig=fig, name="glène")
fig = pltm.add_points(points=p_coracoide_stl, fig=fig, name="coracoide")
fig = pltm.add_points(points=p_acromion_stl, fig=fig, name="acromion")
fig.show()

The same points are easily found on the model

1.1.3. Computations¶

1.1.3.1. Get scanned points from our system scan device¶

df_scan = pd.read_pickle("data_scan.p")
df_scan

	p_glene_scapula			p_acromion_scapula			p_coracoide_scapula			p_scan_scapula
	p3d			p3d			p3d			p3d
	x	y	z	x	y	z	x	y	z	x	y	z
0	0.014046	0.010801	0.118422	-0.023197	-0.022039	0.086035	-0.002085	-0.037214	0.123242	0.015003	-0.004692	0.11943
1	0.014152	0.01073	0.118276	-0.023097	-0.022086	0.085978	-0.002091	-0.037263	0.123189	0.015018	-0.004673	0.119397
2	0.014185	0.010854	0.11838	-0.023238	-0.022144	0.086024	-0.002113	-0.037188	0.123366	0.015072	-0.004562	0.119343
3	0.014362	0.01094	0.118554	-0.023217	-0.022149	0.086317	-0.001946	-0.037163	0.123109	0.015048	-0.004463	0.119347
4	0.014379	0.010899	0.118494	-0.023224	-0.022023	0.086214	-0.001959	-0.036706	0.12319	0.015013	-0.004507	0.11929
...	...	...	...	...	...	...	...	...	...	...	...	...
79	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.015566	-0.001602	0.126646
80	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.01562	-0.001514	0.126536
81	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.015541	-0.001624	0.126567
82	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.015546	-0.001497	0.126532
83	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.015578	-0.001497	0.126524

84 rows × 12 columns

1.1.3.2. Plot scanned points in their reference frame¶

Here the points are expressed in scapula reference frame

cols = set([cols[0] for cols in df_scan.columns])
fig = go.Figure()
for i, c in enumerate(cols):
    fig.add_trace(
        go.Scatter3d(
            x=df_scan[c].p3d.x,
            y=df_scan[c].p3d.y,
            z=df_scan[c].p3d.z,
            name=c,
            marker_size=2,
            mode="markers",
        )
    )
fig.show()

1.1.4. Create initial guess¶

1.1.4.1. Gathering input data of reference points¶

p_glene_scan = np.array(df_scan["p_glene_scapula"].p3d.mean())
p_coracoide_scan = np.array(df_scan["p_coracoide_scapula"].p3d.mean())
p_acromion_scan = np.array(df_scan["p_acromion_scapula"].p3d.mean())

tri_point_stl = np.concatenate([p_glene_stl, p_coracoide_stl, p_acromion_stl]).reshape(
    -1, 3
)
tri_point_scan = np.concatenate(
    [p_glene_scan, p_coracoide_scan, p_acromion_scan]
).reshape(-1, 3)

1.1.4.2. Compute initial guess¶

def unpack(X):
    return X.reshape(2, 3)


def cost(X):
    Rsc2stl, Tsc2stl = unpack(X)
    tri_point_out = gbu.utils.apply_RT(tri_point_scan, Rsc2stl, Tsc2stl)
    dist = abs(tri_point_stl - tri_point_out)
    pen_pin_radius = 0
    res = dist - pen_pin_radius

    return res.flatten()


X0 = np.zeros(6)
sol = optimize.least_squares(
    cost,
    X0,
    method="lm",
    ftol=1.0e-12,
    xtol=1.0e-12,
    gtol=1.0e-10,
)

X_pre = sol.x
R_inital_guess, T_initial_guess = unpack(X_pre)
print(f"Initial guess solution: rvec={R_inital_guess}, tvec={T_initial_guess}")

Initial guess solution: rvec=[ 2.23440972 -0.14576349  0.55625112], tvec=[-0.02428643  0.12368749  0.09225085]

LIFEWORK: coLlection of thEsis WORKs

Surface Matching : example on human scapula

Contents